'''
Utilities to read RINEX GPS file, supports Rinex and Crinex files 
(*.Z *.gz *.tar.gz).

Mostly you will use read_rinex(file) function, a GPSData object is returned.

Based in gpsdata module created by Nick Matteo <kundor@kundor.org>

@author: Juan C. Espinoza
@contact: jucar.espinoza@gmail.com

December 09, 2009

CHANGELOG:

version 1.7 (28/10/2013)
- Fix several bug in merge and save methods od GPSData for supporting GPS L5
  frequency.

version 1.6 (30/07/2013)
- Add TECData and TECrecord classes used to parse TEC files.

version 1.5 (30/06/2012)
- Add merge method to GPSData object.
- Add getgaps method to GPSData object.
- Add plot and plot_tec method to GPSData object.

version 1.2
- Added rinex conversion functions from old gpsProcess module
- Add calcpos method to GPSData object (calculate receiver position)

version 1
- Initial release
'''
import Gps

import os, re, tarfile, itertools
import numpy as np
from cStringIO import StringIO
from subprocess import Popen, PIPE
from math import sin, cos, asin, acos

import Utility as utl
from Gps import *
from headers import *
from Plotter import plot_tec, plot_gpsdata

__version__ = '1.7'

my_float = utl.my_float
my_str   = utl.my_str

rnx_codes = ('C1','C2','P1','P2','L1','L2','S1','S2')
tec_vars = ['eqTEC', 'ltime', 'ele', 'lat', 'lon', 'sTEC', 'azi']

mystrip   = lambda s  : s.strip()
truth     = lambda x  : True
to3float  = lambda s  : tuple([tofloat(s[k*14:(k+1)*14]) for k in (0,1,2)])


def btog(c):
    if c in (None, '', ' '):
        return 'G'
    return c.upper()

def toint(x):
    if x is None or x.strip() == '':
        return 0
    return int(x)

def tofloat(x):
    if x is None or x.strip() == '':
        return 0.
    return float(x)

def choose(a, b):
    if a is not None and b in (' ', None):
        return a
    return b.replace('&', ' ')

def value(thing):
    '''
    Ensure that arbitrary attributes can be set on `thing'.
    E.g. foo = value(foo); foo.bar = 'qux'
    '''

    if type(thing)in (float, int):
        return my_float(thing)
    elif type(thing)==str:
        thing = my_str(thing)
    return thing

class MkRinex(object):
    '''
    Class useful to create a rinex file
    '''

    def __init__(self, marker, interval, obscodes):
        self.marker   = marker
        self.interval = interval
        self.obscodes = obscodes

    def mk_header(self, firsttime, **kwargs):
        '''
        Return a string with a rinex header

        first_obs = date for the first observation [datetime]
        '''

        program   = '%s.py v%s' % ('lisnUtils', __version__)
        firsttime = firsttime.timetuple()[:6]
        obscodes  = '%6i' % len(self.obscodes)        
        for n,obs in enumerate(self.obscodes):
            if n>1 and n%9==0:
                obscodes += '# / TYPES OF OBSERV\n'
                obscodes += '%6s' % ''
            obscodes += '%6s' % obs        
        
        fmt = '%%%is' % ((9-n%9-1)*6+19) 
        obscodes +=  fmt % '# / TYPES OF OBSERV'
        run_date   = GPSDateTime.now().strftime('%Y/%m/%d %H:%M:%S')

        run_by       = kwargs.get('run_by', lisn_header['run_by'])
        observer     = kwargs.get('observer', lisn_header['observer'])
        agency       = kwargs.get('agency', lisn_header['agency'])
        receivernum  = kwargs.get('receivernum', '')
        receivertype = kwargs.get('receivertype', '')
        receiverver  = kwargs.get('receiverver', '')
        antennanum   = kwargs.get('antennanum', '')
        antennatype  = kwargs.get('antennatype', '')
        xyz          = kwargs.get('xyz', (0, 0, 0))+('APPROX POSITION XYZ', )
        delta        = kwargs.get('delta', (0, 0, 0))+('ANTENNA: DELTA H/E/N', )
        interval     = kwargs.get('interval', self.interval)
        endtime      = kwargs.get('endtime')
        comments     = kwargs.get('comment', '')
        comments     = lisn_comment+comments
        comment      = ''
        for line in StringIO(comments):
            if line:
                comment += '%-60sCOMMENT\n' % line[:60].strip()

        header  = '%9s%38s%33s\n' % ('2.10', 'OBSERVATION DATA    G (GPS)',
                                     'RINEX VERSION / TYPE')
        header += '%-20s%-20s%-20sPGM / RUN BY / DATE\n' % (program[:20],
                                                            run_by[:20],
                                                            run_date)
        header += '%s' % comment
        header += '%-60sMARKER NAME\n' % self.marker
        header += '%-20s%-40sOBSERVER / AGENCY\n' % (observer[:20], agency[:40])
        header += '%-20s%-20s%-20sREC # / TYPE / VERS\n' % (receivernum,
                                                            receivertype,
                                                            receiverver)
        header += '%-20s%-20s%-20sANT # / TYPE\n' % (antennanum, antennatype, '')
        header += '%14.4f%14.4f%14.4f%37s\n' % xyz
        header += '%14.4f%14.4f%14.4f%38s\n' % delta
        header += '%6i%6i%68s\n' % (1, 1, 'WAVELENGTH FACT L1/2')
        header += '%s\n' % obscodes
        header += '%10.3f%58s\n' % (interval, 'INTERVAL')
        header += '%6d%6d%6d%6d%6d%13.7f     GPS         TIME OF FIRST OBS\n' % firsttime
        if endtime:
            endtime = endtime.timetuple()[:6]
            header += '%6d%6d%6d%6d%6d%13.7f     GPS         TIME OF LAST OBS\n' % endtime
        header += '%73s\n' % 'END OF HEADER'
        return header

    def mk_satline(self, epoch, prn_list):
        '''
        Create rinex observation record
        '''

        sat_line = ' %02d %2d %2d %2d %2d%11.7f  0%3d' % (epoch.year%100,
                        epoch.month, epoch.day, epoch.hour, epoch.minute,
                        epoch.second, len(prn_list))
        for i, prn in enumerate(prn_list):
            if i>0 and i%12==0:
                sat_line += '\n % 34s' % prn
            else:
                sat_line += '%s' % prn
        return sat_line + '\n'

    def mk_recordline(self, obsvalues, LLI, SNR):
        '''
        Create a rinex observations given obs value, LLI (lost of lock
        indicator) and SNR flag

        output: rinex record in text format
        '''
        format = ''
        nrec   = 0

        for code in self.obscodes:
            #if code not in rnx_codes:
            #    continue
            if nrec%5==0 and nrec>1:
                format+='\n'
            format += '%14.3f%s%s'
            nrec+=1
        format+='\n'
        values = [x[i] for x in zip(obsvalues, LLI, SNR) for i in range(3)]
        return format % tuple(values)

    def snr_flag(self, rec, code=None):
        '''
        '''

        if code=='L1' and 'S1' in self.obscodes and 'S1' in rec:
            snr = rec['S1']
        elif code=='L2' and 'S2' in self.obscodes and 'S2' in rec:
            snr = rec['S2']
        else:
            return 0
        if snr<0:
            return 0
        flags = [0, 26, 28, 32, 36, 39, 42, 45, 49]
        value = [x for x in flags if snr>=x][-1]
        return flags.index(value)+1

class Hatanaka(object):
    '''
    Class to handled hatanaka compress and uncompress files
    '''
    def __init__(self):
        '''
        '''

        ###check for rinex to crinex command
        if os.path.exists('/usr/local/bin/rnx2crx'):
            self.rnx2crx = '/usr/local/bin/rnx2crx'
        else:
            self.rnx2crx = None

        if os.path.exists('/usr/local/bin/crx2rnx'):
            self.crx2rnx = '/usr/local/bin/crx2rnx'
        else:
            self.crx2rnx = None

        if os.path.exists('/usr/bin/compress'):
            self.zip_cmd = '/usr/bin/compress'
        elif os.path.exists('/usr/local/bin/compress'):
            self.zip_cmd = '/usr/local/bin/compress'
        else:
            self.zip_cmd = None
            print 'compress command not found, please install ncompress '

        if os.path.exists('/usr/bin/uncompress'):
            self.unzip_cmd = '/usr/bin/uncompress'
        elif os.path.exists('/usr/local/bin/uncompress'):
            self.unzip_cmd = '/usr/local/bin/uncompress'
        elif os.path.exists('/bin/uncompress'):
            self.unzip_cmd = '/bin/uncompress'
        else:
            self.unzip_cmd = None
            print 'uncompress command not found, please install ncompress '

    def compress(self, file, Z=False, rem=False):
        '''
        Create a CRINEX file from RINEX
        '''
        if self.rnx2crx:
            print 'Creating crinex file: %s' % file[:-1]+'d'
            os.system('cat %s | %s - > %s' % (file, self.rnx2crx,
                                                file[:-1]+'d'))
            if os.path.exists(file[:-1]+'d'):
                new_file = file[:-1]+'d'
                if Z and self.zip_cmd:
                    print 'Compressing file : %s.Z' % new_file
                    os.system('%s -f %s' % (self.zip_cmd, new_file))
                    new_file += '.Z'
                if rem:
                    os.remove(file)
                return new_file
            else:
                return file
        else:
            print 'RNX2CRX not found'
            return file

    def uncompress(self, file):
        '''
        Create RINEX file from CRINEX
        '''
        
        if file.endswith('.Z') and self.unzip_cmd:
            os.system('%s -f %s' % (self.unzip_cmd, file))
            if os.path.exists(file[:-2]):
                file = file[:-2]
                
        if self.crx2rnx and os.path.exists(file) and file.endswith('d'):
            print 'Creating rinex file: %s' % file[:-1]+'o'
            os.system('cat %s | %s - > %s' % (file, self.crx2rnx,
                                                file[:-1]+'o'))
            if os.path.exists(file[:-1]+'o'):
                return file
        else:
            print 'Invalid file: %s' % file
            return False

htnk = Hatanaka()

def versioncheck(ver):
    '''
    Given RINEX format version ver, verify that this program can handle it.
    '''
    nums = ver.split('.')
    if not 0 < len(nums) < 3:
        raise ValueError('RINEX Version not parsable')
    elif int(nums[0]) != 2:
        raise IOError('RINEX File not version 2; unsupported')
    elif len(nums) > 1 and int(nums[1]) > 11:
        warn('RINEX minor version more recent than program.')
    return ver.strip()


def crxcheck(ver):
    '''
    Check whether Compact RINEX version is known to this program.
    '''
    if ver != '1.0':
        raise ValueError('CRINEX version ' + ver + ' not supported.')
    return ver.strip()


def iso(c):
    '''
    Ensure that the character c is `O' (for RINEX observation data.)
    '''
    if c.upper() != 'O':
        raise IOError('RINEX File is not observation data')
    return c.upper()

def headertime(s):
    '''
    Parse RINEX header times into GPSDateTime objects
    '''
    if not s.strip():
        return None
    year   = int(s[0:6])
    if year<1900: year += 1900
    month  = toint(s[6:12])
    day    = toint(s[12:18])
    hour   = toint(s[18:24])
    minute = toint(s[24:30])
    second = tofloat(s[30:43])
    usec   = (second - int(second)) * 1000000
    return GPSDateTime(year, month, day, hour, minute, int(second), int(usec))

def recordtime(s, baseyear=None):
    '''
    Parse RINEX time epoch into GPSDateTime objects
    the latter has two digit years which can be disambiguation with `baseyear'.
    '''
    
    if not s.strip():
        return None
    year = int(s[0:3])
    if baseyear is not None:
        year += (int(baseyear)/100)*100
    elif year < 80:
        year += 2000
    else:
        year += 1900
    month  = toint(s[3:6])
    day    = toint(s[6:9])
    hour   = toint(s[9:12])
    minute = toint(s[12:15])
    second = tofloat(s[15:26])
    usec   = (second - int(second)) * 1000000
    return GPSDateTime(year, month, day, hour, minute, int(second), int(usec))

class wavelength(object):
    '''
    Parse RINEX WAVELENGTH FACT L1/2 headers

    These headers specify 1: Full cycle ambiguities (default),
    2: half cycle ambiguities (squaring), or 0: does not apply,
    either globally or for particular satellites.
    This is only valid for GPS satellites on frequencies L1 or L2.
    '''
    # If prn list is empty (numsats = 0), L1/2 ambiguity applies to all
    # satellites.  Otherwise, it applies to satellites given in the prnlist;
    # continuation lines are allowed.
    # Ambiguity information is valid until the next 'global' WAVELENGTH FACT,
    # or until that prn is reset.
    def __init__(self):
        '''Set all satellites to default wavelength ambiguity, 1.'''
        self.waveinfo = dict([('G%02d' % prn, (1, 1)) for prn in range(1, 33)])

    def __call__(self, s):
        '''Update wavelength ambiguities with information from a new header.'''
        l1amb   = toint(s[0:6])
        l2amb   = toint(s[6:12])
        numsats = toint(s[12:18])
        if not numsats:  # This is a `global' line
            self.waveinfo = dict([('G%02d' % prn, (l1amb, l2amb))
                for prn in range(1, 33)])
        else:
            for p in range(numsats):
                prn = btog(s[21 + 6 * p]) + '%02d' % \
                                              toint(s[22 + 6 * p : 24 + 6 * p])
                self.waveinfo[prn] = (l1amb, l2amb)
        return self.waveinfo.copy()


class obscode(object):
    '''
    Parse RINEX # / TYPES OF OBSERV headers, specifying observation types.

    These header list observation codes which will be listed in this file.
    Continuation lines are necessary for more than 9 observation types.
    It is possible to redefine this list in the course of a file.
    '''
    # There must be `numtypes' many observation codes, possibly over two lines.
    # Continuation lines have blank `numtypes'.
    def __init__(self):
        self.numtypes = None

    def __call__(self, s):
        nt = toint(s[0:6])
        if self.numtypes is not None and not nt:  # continuation line
            if len(self.obstypes) >= self.numtypes:
                raise RuntimeError('Observation code headers seem broken.')
            for ot in range(min(self.numtypes - len(self.obstypes), 9)):
                self.obstypes += [s[6 * ot + 10 : 6 * ot + 12]]
        elif nt:
            self.numtypes = nt
            self.obstypes = []
            for ot in range(min(nt, 9)):
                self.obstypes += [s[6 * ot + 10 : 6 * ot + 12]]
        else:
            raise RuntimeError('Observation type code continuation header '
                               'without beginning!')
        return self.obstypes[:]


class satnumobs(object):
    '''
    Parse RINEX PRN / # OF OBS headers.

    These headers list how many of each observation type were recorded for
    each satellite included in the file.  If present, there will be one for
    each satellite in the file (as reported in the # OF SATELLITES header.)
    If there are more than 9 observation types, a continuation line will be
    necessary for each satellite.
    This program will determine this information anyway, and check against the
    header if it is supplied.
    '''
    def __init__(self):
        self.sno = {}
        self.prn = None

    def __call__(self, s):
        '''Return a dictionary, by satellite PRN code, of observation counts.

        The counts are a list in the same order as obscode().
        '''
        prn = s[0:3]
        if prn.strip() == '' and self.prn is not None:  # continuation line
            pass
        elif prn.strip() != '':
            prn = btog(prn[0]) + '%02d' % toint(prn[1:])
            self.prn = prn
            if prn in self.sno:
                warn('Repeated # OF OBS for PRN ' + prn + ', why?')
            else:
                self.sno[prn] = []
        else:
            raise RuntimeError('PRN / # OF OBS continuation without beginning!')
        for no in range(9):
            obs = s[no * 6 + 3: no * 6 + 9]
            if obs.strip() == '':
                break
            else:
                self.sno[self.prn] += [toint(obs)]
        return self.sno

class SatBias(dict):
    '''
    Read (download if necesary) satellite biases from UNIBE and create a
    dictionary by prn, bias values.
    '''

    ftp_add = 'ftp.unibe.ch'

    def __init__(self, date, obscodes=['C1'], path=None, factor=-TECUns):

        self.date = date
        if not path:
            path = os.path.join(localpath, 'biases')
        if not os.path.exists(path):
            os.makedirs(path)
        self.p1c1 = os.path.join(path, date.strftime('p1c1%y%m.dcb'))
        self.p1p2 = os.path.join(path, date.strftime('p1p2%y%m.dcb'))

        try:
            f_p1p2 = open(self.p1p2)
            f_p1c1 = open(self.p1c1)
        except IOError:
            print 'Attempt to download DCB files...'
            if self.download():
                f_p1p2 = open(self.p1p2)
                f_p1c1 = open(self.p1c1)
            else:
                #raise RuntimeError('Satellite bias files not available.')
                print 'Satellite bias files not available.'
                self.default()
                return

        print 'Satellite Bias files:'
        print self.p1p2
        print self.p1c1

        for line in f_p1p2:
            if line[0] == 'G':
                self[line[:3]] = float([x for x in line.split(' ') \
                                        if len(x)<>0][1])*factor

        if 'C1' in obscodes:
            for line in f_p1c1:
                if line[0] == 'G':
                    self[line[:3]] -= float([x for x in line.split(' ') \
                                             if len(x)<>0][1])*factor

    def download(self):
        '''
        Try to donwload bias files from unibe
        '''
        ftp_path = '/aiub/CODE/'
        actual = datetime.now().year == self.date.year and \
                datetime.now().month == self.date.month
        if not actual:
            ftp_path = os.path.join(ftp_path, `self.date.year`)
            label = '%02d%02d' % (self.date.year%100, self.date.month)
            ext = '.Z'
        else:
            label = ''
            ext = ''

        try:
            ftp = utl.FTP(self.ftp_add)
            ftp.login()
            ftp.cwd(ftp_path)
            ftp.retrbinary('RETR P1C1%s.DCB%s' % (label, ext),
                               open('%s%s' % (self.p1c1, ext), 'wb').write)
            ftp.retrbinary('RETR P1P2%s.DCB%s' % (label, ext),
                               open('%s%s' % (self.p1p2, ext), 'wb').write)
            ftp.close()
        except:            
            print 'Error downloading files'
            return False
        if not actual:
            #raise RuntimeError('Extract downloaded files %s.Z, %s.Z and try again' % \
            #        (self.p1p2, self.p1c1))
            #buf = StringIO(subprocess.check_output(['zcat', filename]))
            os.system('gunzip %s%s' % (self.p1p2, ext))
            os.system('gunzip %s%s' % (self.p1c1, ext))
        return True

    def default(self):
        '''
        '''
        for i in range(33):
            self['G%02d' % i] = 0

class GPSData(list):
    '''
    A GPSData object is primarily a list of records, one for each epoch.
    in chronological order; each record is a dictionary by satellite id
    of dictionaries by observation code of values.

    GPSData.header['name'] also gives access to RINEX header values.
    '''
    def __init__(self):
        self.header    = {}
        self.satsystem = None
        self.date      = None
        self.dates     = {}
        self.rec_bias  = 0
        self.tzoffset  = 0
        self.prns      = set()
        self.inmotion  = False
        self.phasearcs = {}
        self.intervals = set()
        self.station   = {}

    def __str__(self):
        '''
        #TODO
        '''
        types = {'O': 'Observation', 'N':'Navigation'}
        s  = 'Rinex %s Data\n' % types[self.header['filetype']]
        s += ' marker    : %s\n' % self.header['marker'][0]
        s += ' obscodes  : %s\n' % self.obscodes()
        s += ' # records : %d\n\n' % len(self)


        prns = self.phasearcs.keys()
        prns.sort()
        for prn in prns:
            s += ' PRN %s:\n' % prn
            for arc in self.phasearcs[prn]:
                s += '  %s - %s\n' % (self[arc[0]].epoch, self[arc[1]].epoch)

        return s

    def update(self, prn, obs):
        '''
        Add new observables info of a prn to last record.
        '''
        self[-1][prn] = obs        

    def obscodes(self, index=-1):
        '''
        Return (current) list of observation codes stored in this GPSData.
        '''
        if 'obscodes' not in self.header:
            raise RuntimeError('RINEX file did not define data records')

        obs = self.header['obscodes']

        if index == 0:
            index = min(obs)
        elif index == -1:
            index = max(obs)
        else:
            index = max([k for k in obs if k <= index])
        return obs[index]

    def endphase(self, prn):
        '''
        End current phase-connected-arc, if any, for satellite prn.
        Ends arc just before the current record.
        '''
        if isinstance(prn, (list, tuple, set, dict)):
            [self.endphase(p) for p in prn]
        elif prn in self.phasearcs and self.phasearcs[prn][-1][1] is None:
            self.phasearcs[prn][-1][1] = len(self)-1

    def breakphase(self, prn):
        '''
        Begin new phase-connected-arc for satellite prn.
        '''
        if isinstance(prn, (list, tuple, set, dict)):
            [self.breakphase(p) for p in prn]
        elif prn not in self.phasearcs:
            self.phasearcs[prn] = [[len(self) - 1, None]]
        else:
            self.endphase(prn)
            self.phasearcs[prn] += [[len(self) - 1, None]]

    def checkbreak(self, slip=False):
        '''
        Check times slip and bad records slips.

        Checks at last record added.  This should be called for each record
        inserted, after all its values have been added.

        Also update attributes dates, prns, intervals
        '''

        #update dates
        date = self[-1].epoch.date
        if date not in self.dates:
            self.dates[date]  = 1
        else:
            self.dates[date] += 1
        #update prns
        self.prns.update(self[-1].keys())
        #update intervals
        if len(self)>1:
            self.intervals.add(self[-1].epoch-self[-2].epoch)
        #check for breaks
        for prn in set(self[-1]).union(self.phasearcs):
            if slip:
                self.breakphase(prn)
                continue

            if len(self)>1:
                if self[-1].epoch-self[-2].epoch>1800:
                    self.breakphase(prn)
                    continue

            bad = self[-1].badness(prn)

            if not bad and prn not in self.phasearcs:
                self.phasearcs[prn] = [[len(self)-1, None]]
                continue
            if not bad and self.phasearcs[prn][-1][1] is not None:
                self.phasearcs[prn] += [[len(self)-1, None]]
                continue
            if prn not in self.phasearcs or self.phasearcs[prn][-1][1] is not \
                    None:
                continue
            if bad:
                self.endphase(prn)

    def timesetup(self, filename=None):
        '''
        Add significant attributes to object, also set up year disambiguation.
        '''
        if self.satsystem is None and 'satsystem' in self.header:
            self.satsystem = self.header['satsystem']
        baseyear = None
        if 'firsttime' in self.header:
            baseyear = self.header['firsttime'].year
        if 'endtime' in self.header:
            if baseyear is None:
                baseyear = self.header['endtime'].year
        '''
        if filename:
            filename = filename.split('/')[-1]
            if re.search('[0-9]{6}', filename):
                fdate = time.strptime(re.search('[0-9]{6}', filename).group(),
                    '%y%m%d')
                self.date = GPSDateTime(*fdate[:3])
        '''
        return baseyear

    def iterlist(self, prn, obscode=None, filter=None):
        '''
        Returns an iterator over the list of records given, prn and obscodes.
        '''
        def chooser(obs, rec, prn):
            if obs == 'epoch':
                return rec['epoch']
            elif obs in rec[prn]:
                return rec[prn][obs]
            else:
                return None

        rec_idx = []
        for arc in self.phasearcs[prn]:
            if filter:
                rec_idx += [x for x in range(arc[0], arc[1]) \
                            if self[x][prn][filter[0]]>filter[1]]
            else:
                rec_idx += [x for x in range(arc[0], arc[1])]

        if isinstance(obscode, (list, tuple, set, dict)):
            if not obscode:
                obscode = None
            elif isinstance(obscode, (tuple, set, dict)):
                obscode = [o for o in obscode]
        elif isinstance(obscode, str):
            obscode = [obscode]
        if obscode is None:
            obscode = list(self.allobs)

        for x in rec_idx:
            yield [chooser(obs, self[x], prn) for obs in obscode]

    def check(self):
        '''
        Check header values of file with found and correct them.
        Also save useful attributes.
        '''        
        #checking headers values
        if 'interval' in self.header:
            if self.header['interval'] != min(self.intervals):
                warn('INTERVAL ' + str(self.header['interval']) + \
                        ' does not match minimum observation interval ' +\
                        str(min(self.intervals)))
                self.header['interval'] = min(self.intervals)
        else:
            self.header['interval'] = min(self.intervals)

        if self.header['interval']<=0:
            warn('Zero or negative interval found...check')
            long = 1800/30
        else:
            long = 1800/self.header['interval']
        for prn in self.phasearcs:
            self.phasearcs[prn] = [arc for arc in self.phasearcs[prn] \
                                   if arc[1]-arc[0]>long]

        #save date attribute
        N = 0
        for dt,n in self.dates.items():
            if n>N:
                N = n
                self.date = dt
        #save tzoffset attribute
        if 'longitude' in self.station:
            self.tzoffset = self.station['longitude']/15.
        else:
            lat, lon, height = xyz2lla(*self.header['markerpos'][0])
            self.tzoffset = lon/15.

    def merge(self, data_i):
        '''
        Merge current gpsdata object with another.
        '''

        data_j = self[:]

        self.dates     = {}
        self.rec_bias  = 0
        self.prns      = set()
        self.phasearcs = {}
        self.intervals = set()
        
        for obs in data_i.header['obscodes'][0]:
            if obs not in self.header['obscodes'][0]:
                self.header['obscodes'][0].append(obs)
        
        del self[:]
        
        i = 0
        j = 0
        imax = len(data_i)
        jmax = len(data_j)

        while True:          
            if data_i[i]==data_j[j]:
                self.append(data_j[j])
                i+=1
                j+=1
            elif data_i[i]>data_j[j]:
                self.append(data_j[j])
                j+=1
            else:
                self.append(data_i[i])
                i+=1
            self.checkbreak()            

            if i==imax or j==jmax:
                break
        if i==imax:
            temp = data_j[j:]
        else:
            temp = data_i[i:]

        for rec in temp:
            self.append(rec)
            self.checkbreak()            
        self.endphase(self.prns)
        self.check()
        print 'Merging old:%d + new:%d -> %d records' % (jmax, imax, len(self))

    def getgaps(self):
        '''
        Return a list of gaps presents
        '''
        interval = min(self.intervals)
        gaps  = []
        for i in range(1,len(self)):
            delta = self[i].epoch-self[i-1].epoch
            if delta>interval:
                gaps.append([self[i-1].epoch, self[i].epoch, delta])
        return gaps

    def calcpos(self, obscode=None, local=False, samples=20):
        '''
        Calculate receiver position using least squares method,
        -- NO CORRECTIONS ARE APPLIED TO SOLUTION --
        '''

        if obscode:
            obscode = obscode.upper()
        else:
            if 'C1' in self.header['obscodes'][0]:
                obscode = 'C1'
            elif 'P2' in self.header['obscodes'][0]:
                obscode = 'P2'
            else:
                obscode = 'P1'
        if obscode not in self.header['obscodes'][0]:
            raise RuntimeError('Obscode "%s" missing: receiver position can not be calculate' % obscode)
        
        ephe = Yuma(self.date, local=local)
        
        ret  = []
        step = len(self)/samples

        if 'location' in self.station:
            xyz = lla2xyz(*self.station['location'])
        else:
            xyz = self.header['markerpos'][0]

        for i in range(0, len(self), step):
            rec = self[i]
            prns = [prn for prn in rec.keys() if prn in ephe]
            if len(prns)<4: continue
            XYZ_rec = np.tile(lla2xyz(*xyz), (len(prns), 1))
            for j in xrange(3):
                #calculate satellite position
                XYZ_sat  = np.array(ephe.sat_pos(prns, rec.epoch, XYZ_rec[0]))
                R_sat    = ((XYZ_sat[:,0]-XYZ_rec[:,0])**2+(XYZ_sat[:,1]-XYZ_rec[:,1])**2+(XYZ_sat[:,2]-XYZ_rec[:,2])**2)**0.5
                R_sat_i  = np.tile(R_sat, (3, 1)).T
                #form matrix system
                A = np.column_stack([(XYZ_sat-XYZ_rec)/R_sat_i, np.ones(len(prns))])
                Y = R_sat-np.array(rec.iterlist(prns, obscode))-np.array(ephe.get_sat_clock_correction(prns, rec.epoch))
                #solve linear system
                DXYZ = np.linalg.solve(np.dot(A.T, A), np.dot(A.T, Y))
                #recalculate reciver position
                XYZ_rec += np.tile(DXYZ[:3], (len(prns), 1))
            ret.append(XYZ_rec[0])
            #break
        ret = np.array(ret)
        #print ret[:,0], ret[:,1], ret[:,2]
        return ret[:,0].mean(), ret[:,1].mean(), ret[:,2].mean()

    def calctec(self, level=5, rec_bias=None, bad_prn=[], path=localpath):
        '''
        Calculate TEC, STEC and eqTEC and append as observation to records.
        the calculation is only for those records in phasearcs.
        level:
            0 = only calculate atec rtec ele azi lat lon
            1 = correct bad points
            2 = fix jumps
            3 = correct cycle slips #TODO
            4 = calc RX bias and eqTEC
            5 = leveling zero TEC
        '''
        if len(self)<50 or not self.phasearcs:
            raise RuntimeError('Not enough data to calculate TEC')
        self.bad_prn = bad_prn
        print 'Calculating TEC...'

        elem = gpselaz(self.date, os.path.join(path, 'almanac'))
        if 'location' in self.station:
            elem.update(*self.station['location'])
        else:
            elem.update(*xyz2lla(*self.header['markerpos'][0]))
        nprns = []
        for prn in self.phasearcs:
            if prn not in elem or not self.phasearcs[prn]:
                nprns.append(prn)
        for prn in nprns:
            del self.phasearcs[prn]

        print 'Calculating: aTEC, rTEC, ele, azi, lat and lon'
        for prn in self.phasearcs.keys():
            for x in self.get_index(prn):
                if prn not in self[x]: continue
                [xazi, xele, xlat, xlon] = elem.find(prn, self[x].epoch.hours)
                self[x][prn]['ele'] = xele
                self[x][prn]['azi'] = xazi
                self[x][prn]['lat'] = xlat
                self[x][prn]['lon'] = xlon
                self[x][prn]['ltime'] = self[x].epoch.hours + xlon/15.
                self[x][prn]['rTEC'] = self[x].ptec(prn)
                self[x][prn]['aTEC'] = self[x].ctec(prn)

        if level>=1:
            print 'Correcting bad points'
            for prn, arclist in self.phasearcs.items():
                for arc in arclist:
                    [self.correct_points(prn, x) for x in range(*arc)[1:-1]]
        if level>=2:
            print 'Fixing jumps'
            for prn, arclist in self.phasearcs.items():
                for arc in arclist:
                    rtec = [self[x][prn]['rTEC'] for x in range(*arc)]
                    times = [self[x].epoch.doys for x in range(*arc)]
                    new_rtec = utl.fix_jumps(rtec, times, 1.173, id=prn)
                    for x in range(len(new_rtec)):
                        self[x+arc[0]][prn]['rTEC'] = new_rtec[x]
        if level>=3:
            print 'Processing data Cycle Slips'
            sat_bias = SatBias(self.date, self.obscodes(0), os.path.join(path,
                                                                    'biases'))
            #for bprn in bad_prn:
            #    delta = self.correct_sat_bias(bprn)
            #    if delta: sat_bias[bprn] += delta
            #sat_bias['G32'] += 10
            #sat_bias['G22'] -= 7
            self.calctec2(self.prns, sat_bias, 0)#, True)
            for prn, arclist in self.phasearcs.items():
                for arc in arclist:
                    self.correct_slips(prn, arc)
        if level>=4:
            bias_file = '%s.bias' % self.header['marker'][0]
            bias_file = os.path.join(path, 'biases', bias_file)
            bias_list = {}
            old_rec_bias = None
            write_bias = True
            leveling = True            
            if os.path.exists(bias_file):                
                for line in open(bias_file):
                    bias_list[GPSDateTime(line.split(',')[0], fmt='%Y %m %d')] = float(line.split(',')[1].strip())            
            if rec_bias:
                print 'Calculating eqTEC using RX bias =', rec_bias
                self.rec_bias = rec_bias
                self.calctec2(self.prns, sat_bias, self.rec_bias)
                if level>=5:
                    print 'Leveling bias for negative TEC'
                    self.zerolevel(sat_bias)
                    if abs(rec_bias-self.rec_bias)>3:
                        print 'zero leveling bias difference to high using=%s' % rec_bias
                        self.rec_bias = rec_bias
                        self.calctec2(self.prns, sat_bias, self.rec_bias)
            else:
                print 'Calculating RX bias'
                self.rec_bias = self.rx_bias(sat_bias)
                
                if bias_list:
                    dt_list = [dt for dt in bias_list]
                    dt_list.sort()
                    for dt_bias in dt_list:
                        if dt_bias<self.date:
                            old_rec_bias = bias_list[dt_bias]
                            old_dt = dt_bias
                        else:
                            break
                if self.rec_bias and old_rec_bias:
                    if abs(old_rec_bias-self.rec_bias)>10:
                        warn('Using RX bias=%.1f, date: %s' % (old_rec_bias,
                                                               dt))
                        self.rec_bias = old_rec_bias
                        #write_bias = False
                elif not self.rec_bias and old_rec_bias:
                    warn('Using RX bias=%.1f, date: %s' % (old_rec_bias, dt))
                    self.rec_bias = old_rec_bias
                    write_bias = False
                    leveling = False
                elif not self.rec_bias and not old_rec_bias:
                    raise RuntimeError('RX bias can not be estimated')
                self.calctec2(self.prns, sat_bias, self.rec_bias)
                if level>=5 and leveling:
                    print 'Leveling bias for negative TEC'
                    self.zerolevel(sat_bias)
                    #if old_rec_bias and abs(old_rec_bias-self.rec_bias)>8:
                    #    self.rec_bias = old_rec_bias

                if write_bias:
                    print 'Updating Receiver Bias file: %s' % bias_file
                    bias_list[self.date] = self.rec_bias
                    dt_list = [dt for dt in bias_list]
                    dt_list.sort()
                    fbias = open(bias_file, 'w')
                    [fbias.write('%d %02d %02d, %.1f\n' % (dt.year,
                        dt.month, dt.day, bias_list[dt])) for dt in dt_list]
                    fbias.close()
            #self.correct_sat_bias('G32')

    def correct_points(self, prn, x):
        '''
        Correct points too distant
        '''
        if prn not in self[x-1] or prn not in self[x] or prn not in self[x+1]:
            return
        if abs(self[x][prn]['rTEC']-self[x-1][prn]['rTEC'])>1.173:
            if abs(self[x+1][prn]['rTEC']-self[x-1][prn]['rTEC'])<1.173:
                print '1 point corrected', prn, self[x].epoch
                self[x][prn]['rTEC'] = 0.5*(self[x+1][prn]['rTEC'] + \
                                        self[x-1][prn]['rTEC'])

    def badrecords(self, prn, ttype, tmax, tdiff):
        '''
        Return a index list with bad values for tec calculation
        '''
        bad = set()
        for arc in self.phasearcs[prn]:
            rec_idx = range(*arc)
            bad.union(set([x for x in rec_idx if self[x][prn][ttype]>tmax]))
            tec_tmp = [self[x][prn][ttype] for x in rec_idx if x not in bad]

            if len(tec_tmp)<>0:
                ave = utl.Stats(tec_tmp).average
            else:
                continue

            bad.union(set([x for x in rec_idx if abs(self[x][prn][ttype]-ave)>tdiff]))
            tmp = [self[x][prn][ttype] for x in rec_idx \
                   if x not in bad]

            if len(tmp)<>0:
                stats = utl.Stats(tmp)
                ave   = stats.average
                stdev = stats.std_dev
                s1    = ave+stdev*4.
                s2    = ave-stdev*4.
                bad.union(set([x for x in rec_idx if (s1<self[x][prn][ttype]) \
                           or (self[x][prn][ttype]<s2)]))
        return bad

    def correct_slips(self, prn, arc):
        '''
        Detect and correct cycle slips taken 10 continues samples and ...
        '''
        return

        nsamp = 6
        rec_idx = range(arc[0]+1, arc[1])
        for x in rec_idx[1:]:
            s1 = x-nsamp/2
            if s1<arc[0]+1:
                s1=arc[0]+1
            s2 = nsamp+s1
            if s2>arc[1]:
                s2 = arc[1]
                s1 = s2-nsamp
            pdiff = [self[x][prn]['rTEC']-self[x-1][prn]['rTEC'] \
                     for x in range(s1,x)+range(x,s2)]
            cdiff = [self[x][prn]['aTEC'] for x in range(s1,x)+range(x,s2)]
            pstat = utl.Stats(pdiff)
            cstat = utl.Stats(cdiff)
            tdiff = self[x][prn]['rTEC']-self[x-1][prn]['rTEC']
            if abs(tdiff)>1.173:
                #print prn, self[x][prn]['rTEC'], tdiff, pdiff, cdiff
                plim_dw = pstat.average-4.*pstat.std_dev
                plim_up = pstat.average+4.*pstat.std_dev
                if plim_up<tdiff or tdiff<plim_dw:
                    #potential cycle slip occurr
                    gdiff = abs(self[x][prn]['TEC'] - self[x][prn]['aTEC'])
                    clim_up = cstat.average + cstat.std_dev
                    clim_dw = cstat.average - cstat.std_dev
                    if clim_dw<self[x][prn]['TEC']<clim_up and abs(tdiff)>1.9:
                        continue
                    if self[x][prn]['ele']<30:
                        #fix as a jump
                        rtec = [self[k][prn]['rTEC'] for k in range(*arc)]
                        secs = [self[k].epoch.seconds for k in range(*arc)]
                        rtec = utl.fix_jumps(rtec, secs, 1.173, x, id=prn)
                        for k in range(x, len(rtec)):
                            self[x+arc[0]][prn]['rTEC'] = rtec[x]
                        continue
                    islips = abs(tdiff/1.173)
                    xadd = 1.173*abs(tdiff)/tdiff*float(islips)
                    print 'Correcting cycleslip prn=%s, idx=%d, tecdif=%.3f' % \
                            (prn, x, tdiff)
                    for y in range(x, arc[1]):
                        self[y][prn]['rTEC'] = self[y][prn]['rTEC'] - xadd

    def ave_sdev(self, prn, obs='eqTEC'):
        '''
        '''
        arcs = self.phasearcs[prn]
        for arc in arcs:
            print prn, arc
            idx_list = range(*arc)
            '''
            prn_list0 = set([self[x].keys()[y] for x in idx_list for y in \
                            range(len(self[x].keys()))])
            print prn_list0
            prn_list = []
            for nprn in prn_list0:
                bad = [self[x].badness(nprn) for x in idx_list]
                if True in bad: continue
                prn_list += [nprn]
            print prn_list
            '''

            eqtec_list = [self[x][nprn][obs] for x in idx_list for nprn \
                          in self[x].keys() if nprn<>prn and obs in self[x][nprn]]
            stat_all = utl.Stats(eqtec_list)
            stat = utl.Stats([self[x][prn][obs] for x in idx_list])
            print stat.average, stat.std_dev
            print stat_all.average, stat_all.std_dev
        return

    def rx_bias(self, sat_bias):
        '''
        Receiver Bias estimation, using mnbrak and brent methods,
        minimizing sum[var('eqTEC')] betwen 03:00 and 06:00 LT.
        '''

        t1 = 3 - self.tzoffset
        t2 = 6 - self.tzoffset

        idx_list = [x for x in range(len(self)) \
                    if t1<self[x].epoch.hours<t2]

        if len(idx_list)==0:
            return False
        else:
            prn_list1 = []
            for x in idx_list:
                prn_list1 += self[x].keys()
            prn_list1 = set(prn_list1)

            prn_list = []
            for prn in prn_list1:
                bad = [self[x].badness(prn) for x in idx_list]
                if True in bad or prn in self.bad_prn or prn not in self.phasearcs:
                    continue
                prn_list += [prn]
            #print prn_list
            #find a range where there is a minimun
            val, fval = utl.mnbrak(100, 50, self.fvarsum,
                                   (idx_list, prn_list, sat_bias))
            #find a minimum of fvarsum function
            bias = utl.brent(val, fval, self.fvarsum,
                             (idx_list, prn_list, sat_bias))
            print 'RX bias = %.2f, sum_var(eqTEC) = %.2f after %d iterations' \
                    % (int(bias[0]), bias[1], bias[2])
        if int(bias[0])==0:
            return 1
        return int(bias[0])

    def correct_sat_bias(self, qprn):
        '''
        '''
        if qprn not in self.phasearcs:
            print 'No correction can be done for prn=%s' % qprn
            return None
        intersc = []
        arcs = self.phasearcs[qprn]
        for arc in arcs:
            #check for an arc bigger than 30 minutes
            if self[arc[1]].epoch-self[arc[0]].epoch>30*60:
                xstr = self[arc[0]].epoch.hours + self.tzoffset
                xend = self[arc[1]].epoch.hours + self.tzoffset + 1
                #self.cls2app()
                xdmin = 5.0
                hdmin = 500.0

                qprn_idx = [x for x in range(*arc) if xstr<self[x][qprn]['ltime']<xend \
                            and self[x][qprn]['ele']>20]

                arc_prn = set([prn for x in qprn_idx for prn in self[x].keys() \
                               if prn<>qprn])

                for aprn in arc_prn:
                    distmn = 10000
                    aprn_idx = [x for x in self.get_index(aprn) if \
                                xstr<self[x][aprn]['ltime']<xend and self[x][aprn]['ele']>20]

                    idxs = [x for x in aprn_idx if x in qprn_idx]

                    for x in idxs:
                        dist = abs(self[x][qprn]['lat']-self[x][aprn]['lat'])
                        if dist < distmn:
                            hdist = self.GTODDBL(self[x][qprn]['lat'], self[x][qprn]['lon'],
                                                 self[x][aprn]['lat'], self[x][aprn]['lon'])
                            distmn = dist
                            idx = x
                    #print qprn, aprn, idx, distmn, hdist
                    if distmn<xdmin and hdist<hdmin:
                        intersc += [(aprn, idx)]
                #####
        #print intersc
        if intersc:
            sum1, sum2, sum3 = 0, 0, 0
            for bub in intersc:
                iprn, x = bub
                sf   = 1.0/cos(asin(0.94092*cos(self[x][qprn]['ele']*deg2rad)))
                sum1 = sum1 + self[x][qprn]['eqTEC']/sf
                sum2 = sum2 + self[x][iprn]['eqTEC']/sf
                sum3 = sum3 + 1.0/sf**2
            bias_crrc = (sum1-sum2)/sum3
            print 'Correcting bias for %s: ' % qprn , bias_crrc
            return bias_crrc
        else:
            print 'No correction can be done for', qprn
            return None

    def get_index(self, prn):
        '''
        Return a list of index given a prn
        '''
        if prn in self.phasearcs:
            return [x for arc in self.phasearcs[prn] for x in range(*arc)]
        else:
            return []

    def GTODDBL(self, XLAT,XLON,SLAT,SLON):
        '''
        '''
        XALT  = 200.0
        STLAT = SLAT*deg2rad
        RE    = 6378.39*(.99832 + .00168*cos(2.*STLAT))
        C     = abs(SLON - XLON)
        AS    = deg2rad*(90.0 - SLAT)
        BS    = deg2rad*(90.0 - XLAT)
        COSC  = cos(AS)*cos(BS) + sin(AS)*sin(BS)*cos(C*deg2rad)
        if COSC >= 0.999999999: COSC = 0.999999999
        if COSC <= -.999999999: COSC = -0.999999999
        HDIST = acos(COSC)*(RE + XALT)
        return abs(HDIST)

    def fvarsum(self, rec_bias, idx_list, prn_list, sat_bias):
        '''
        calculate sum(std_dev(eqTEC)) per sat for a range of epochs 'idx_list'
        '''
        self.calctec2(prn_list, sat_bias, rec_bias)
        sd = 0

        for x in idx_list:
            prns = self[x].keys()
            eqTEC_list = [self[x][prn]['eqTEC'] for prn in prns \
                         if prn in prn_list and 'eqTEC' in self[x][prn]]
            if eqTEC_list:
                sd += utl.Stats(eqTEC_list).std_dev
        return sd


    def calctec2(self, prns, sat_bias, rec_bias):
        '''
        TEC from carrier phase (rTEC) is smooth but ambiguous.
        TEC from code (aTEC), (C1, C2) or encrypted code (P1, P2) is
        absolute but noisy.  Phase data needs to be fitted to pseudorange data.
        This function calculates the TEC for each PRN in each record, where
        possible.
        Also calculate equivalent TEC using satellite bias and receiver bias
        corrections.
        '''
        prns = [prn for prn in prns if prn in self.phasearcs]
        for prn in prns:
            #bads = self.badrecords(prn, 'aTEC', 999, 120)
            #bads.union(self.badrecords(prn, 'rTEC', 9999999, 1200))
            #print prn, bads
            #bads = []
            for arc in self.phasearcs[prn]:
                #suma = [(self[x][prn]['aTEC'] - self[x][prn]['rTEC']) for x in range(arc[0], arc[1]) if x not in bads]
                suma = [(self[x][prn]['aTEC'] - \
                         self[x][prn]['rTEC'])*sin(self[x][prn]['ele']*deg2rad) \
                            for x in range(*arc)]
                xcount = [sin(self[x][prn]['ele']*deg2rad) \
                          for x in range(*arc)]
                avediff = sum(suma)/sum(xcount)
                #avediff = sum(suma)/len(suma)
                for x in range(*arc):
                    self[x][prn]['sTEC'] = self[x][prn]['rTEC'] + avediff
                    self[x][prn]['TEC'] = self[x][prn]['sTEC'] - \
                                                sat_bias[prn]-rec_bias
                    sf = 1.0/cos(asin(0.94092*cos(self[x][prn]['ele']*deg2rad)))
                    self[x][prn]['eqTEC'] = self[x][prn]['TEC']/sf

    def zerolevel(self, sat_bias):
        '''
        Correct receiver bias by leveling zero TEC values.
        '''
        cnt_zero = 0
        cnt_high = 0
        while True:
            prnszero = set()
            for prn in self.phasearcs.keys():
                if prn in self.bad_prn: continue
                try:
                    min_tec = min([self[x][prn]['eqTEC'] for x in self.get_index(prn) if self[x][prn]['ele']>15])
                    if min_tec<0: prnszero.add(prn)
                except ValueError:
                    continue
            if len(prnszero)==0 and cnt_zero==0:
                if cnt_high>20:
                    print 'Max iteration reach'
                    break
                self.rec_bias += 1
                print 'High TEC values found'
                print 'Trying RX bias = %.2f' % self.rec_bias
                self.calctec2(self.prns, sat_bias, self.rec_bias)
                cnt_zero = 0
                cnt_high += 2
            elif len(prnszero)==0:
                break
            elif len(prnszero)>10:
                self.rec_bias -= 5
                print 'Very low TEC'
                print 'Trying RX bias = %.2f' % self.rec_bias
                self.calctec2(self.prns, sat_bias, self.rec_bias)
                cnt_zero = 0
            elif cnt_zero > 6:
                print 'Ending Zero leveling after 7 iterations'
                print 'prns with TEC<0: ', [prn for prn in prnszero]
                break
            else:
                self.rec_bias -= 1
                cnt_zero += 1
                print 'prns with TEC<0: ', [prn for prn in prnszero]
                print 'Trying RX bias = %.2f' % self.rec_bias
                self.calctec2(self.prns, sat_bias, self.rec_bias)

    def vtec(self, max_ele=25):
        '''
        Experimental total vertical tec calculation
        '''
        def calcvTEC(rec, max_ele):
            tec = 0
            count = 0
            for prn in rec.keys():
                if 'eqTEC' in rec[prn] and rec[prn]['ele']>=max_ele and \
                                                prn not in self.bad_prn:
                    tec += rec[prn]['eqTEC']
                    count += 1
            if count:
                return tec/count
            else:
                return 0
        tim = [rec.epoch for rec in self]
        VTEC = [calcvTEC(rec, max_ele) for rec in self]
        #eqTEC = utl.fix_jumps(eqTEC, tim, 3, id='eqTEC')
        return [(tim[x], VTEC[x]) for x in range(len(VTEC)) \
                    if not tim[x].minute%5 and tim[x].second==0]

    def save(self, path='./', marker='code', alt_name=False, hours=None,
             seconds=None, **kwargs):
        '''
        Save data into a new Rinex file use kwargs for header values.
        '''

        def attr(value, atr):
            if hasattr(value, atr): return getattr(value, atr)
            else: return 0
        if hours:
            date = self[-1].epoch - timedelta(1./24*hours)
            if date<self[0].epoch:
                date = self[0].epoch
        else:
            date = self.date

        if alt_name:
            rnx_name = '%s.%02do' % (marker, date.year%100)
        else:
            rnx_name = '%s%003d0.%02do' % (marker, date.doy, date.year%100)

        obscodes = self.obscodes(0)
        rnx = MkRinex(marker, min(self.intervals), obscodes)
        file_tmp = os.tmpfile()
        nrec = 0
        intervals = set()
        for rec in self:
            if rec.epoch<date:
                continue
            if not hours and rec.epoch.day<>date.day:
                break
            if seconds and rec.epoch.second not in seconds:
                continue
            if nrec>0:
                intervals.add(rec.epoch-last_rec.epoch)
            else:
                first_rec = rec
            last_rec = rec
            prns = rec.keys()
            prns.sort()
            sat_line = rnx.mk_satline(rec.epoch, prns)
            dat_line = ''            
            for prn in prns:
                obs_values = [rec[prn][code] for code in obscodes]
                lli_values = [attr(rec[prn][code], 'lli') for code in obscodes]
                snr_values = [attr(rec[prn][code], 'str') for code in obscodes]
                dat_line += rnx.mk_recordline(obs_values, lli_values,
                                snr_values)
            file_tmp.write(sat_line+dat_line)
            endtime = rec.epoch
            nrec += 1

        if nrec<>0:
            kwargs['endtime'] = endtime
            kwargs['interval'] = min(intervals)
            header = rnx.mk_header(first_rec.epoch, **kwargs)
            if not os.path.exists(path):
                os.makedirs(path)
            file_tmp.seek(0,0)
            file_rnx = open(os.path.join(path, rnx_name), 'w')
            file_rnx.write(header)
            file_rnx.write(file_tmp.read())
            file_tmp.close()
            file_rnx.close()
            print 'Rinex records saved %d of %d' % (nrec, len(self))
            return os.path.join(path, rnx_name), nrec
        else:
            file_tmp.close()
            print 'No records found for given date'
            return None, 0

    def save_data(self, vars=tec_vars, path=localpath, prns=None, gz=False):
        '''
        Save observables data into a new column text file.
        '''
        if not os.path.exists(path):
            os.makedirs(path)

        if 'lcode' in self.station:
            lcode = self.station['lcode']
            fullname = self.station['fullname']
            location = self.station['location']
            receiver = self.station['receiver']
        else:
            lcode = self.header['marker'][0]
            fullname = self.header['marker'][0]
            location = Gps.xyz2lla(*self.header['markerpos'][0])
            receiver = self.header['receivertype']

        fname = '%s_%02d%02d%02d.dat' % (lcode, self.date.year%100,
                                         self.date.month, self.date.day)

        fname = os.path.join(path, fname)

        if gz:
            fd = gzip.open(fname+'.gz', 'wb')
            print 'Creating TEC file: %s.gz' % fname
        else:
            fd = open(fname, 'w')
            print 'Creating TEC file: %s' % fname

        header = '%9s%s\n' % (' ', fullname)
        header += '%13.4f%13.4f%13.4f\n' % location
        header += '%s\n' % receiver
        header += 'RINEX\n'
        header += '%5d%5d%5d' % (self.date.year,
                                       self.date.month, self.date.day)
        fd.write(header)
        if prns:
            prns = [prn for prn in prns if prn in self.phasearcs]
        else:
            prns = [prn for prn in self.phasearcs.keys() if prn not in self.bad_prn]
        prns.sort()
        format = '\n%5d    0%4d%4d%4d' + ''.join(['%9.3f' for x in vars])
        for prn in prns:
            nrec = 0
            data = ''
            arcs = self.phasearcs[prn]
            for arc in arcs:
                for x in range(arc[0], arc[1]):
                    #values = [self[x][prn][var] for var in vars]
                    values = []
                    for var in vars:
                        if var in self[x][prn]:
                            values.append(self[x][prn][var])
                        elif var == 'time':
                            values.append(self[x].epoch.hours)
                        else:
                            values.append(-1)
                    nrec += 1
                    data += format % ((nrec, self[x].epoch.hour,
                           self[x].epoch.minute, self[x].epoch.second)\
                           +tuple(values))
            prn_line = '\nPRN %7d%7d' % (int(prn[1:]), nrec)
            fd.write(prn_line+data)
        fd.close()

    def plot_tec(self, figname=None, localtime=False, min_ele=25):
        '''
        Plot TEC
        '''
        
        for prn in self[0]:
            if 'eqTEC' in self[0][prn]:
                return plot_tec(self, figname, localtime)
        warn('TEC not available, run first "caltec" method')

    def plot(self, figname=None, localtime=False):
        '''
        PLOT data availavility
        '''

        plot_gpsdata(self, figname, localtime)

class GPSrecord(dict):
    '''
    A record of observations (many satellites, many channels) at a given epoch.

    Has atributtes epoch, flag and clock offset in addition to a dictionary
    (by PRN code) of dictionaries (by RINEX observation code, e.g. C1, L2)
    of values.
    Can access as record.epoch, record[13], record['G17'], or iteration.
    '''

    def __init__(self, rec, clkoffset):
        self.epoch = rec.epoch
        self.flag = rec.flag
        self.clockoffset = clkoffset

    def __getitem__(self, index):
        '''
        Allow you to access GPS satellites, eg record['G13'], as
        simply record[13].  For GLONASS or Galileo, you must use the full code.
        '''
        if index == 'epoch':
            return self.epoch
        elif isinstance(index, (int, long, float)):
            return dict.__getitem__(self, 'G%02d' % index)
        else:
            return dict.__getitem__(self, index)

    def __contains__(self, index):
        '''
        Allow containment tests (eg if 13 in record:) for abbreviated GPS PRNs.
        '''
        if isinstance(index, (int, long, float)):
            return dict.__contains__(self, 'G%02d' % index)
        return dict.__contains__(self, index)

    def __eq__(self, other):
        '''
        Equal comparisson reference to the epoch of the record
        '''
        
        return self.epoch==other.epoch

    def __ne__(self, other):
        '''
        Not equal comparisson reference to the epoch of the record
        '''

        return self.epoch!=other.epoch

    def __lt__(self, other):
        '''
        Comparisson reference to the epoch of the record
        '''

        return self.epoch<other.epoch

    def __le__(self, other):
        '''
        Comparisson reference to the epoch of the record
        '''

        return self.epoch<=other.epoch

    def __gt__(self, other):
        '''
        Comparisson reference to the epoch of the record
        '''

        return self.epoch>other.epoch

    def __ge__(self, other):
        '''
        Comparisson reference to the epoch of the record
        '''

        return self.epoch>=other.epoch


    def iterlist(self, prns, obscode):
        '''
        '''
        return [self[prn][obscode] for prn in prns]

    def ptec(self, prn):
        '''
        Phase (carrier) TEC, if observations are available.

        Convert cycles of L1, L2 to ns (divide by frequency)
        Then TEC_P (TECU) = (L1(ns) - L2(ns)) * TECU/ns
        1 TECU = 10*16 el / m**2
        Suffers from integer ambiguity (`constant' offset which changes after
        each cycle slip.)  Should be dealt with in `arcs' between cycle slips.
        '''
        L1ns = self[prn]['L1']/F1
        L2ns = self[prn]['L2']/F2
        return (L1ns - L2ns) * TECUns

    def ctec(self, prn):
        '''
        Code (pseudorange) TEC, if observations are available.

        p(ns) = p(m) / .3; TEC_C (TECU) = (p2(ns) - p1(ns) - HWCAL) * TECU/ns
        HWCAL is the hardware calibration in ns
        Suffers from satellite bias, receiver bias, multipath and noise.
        '''
        if 'C1' in self[prn] and 'P2' in self[prn]:
            return (self[prn]['P2'] - self[prn]['C1']) * TECUns/(LIGHT_SPEED*1e-9)
        if 'P1' in self[prn] and 'P2' in self[prn]:
            return (self[prn]['P2'] - self[prn]['P1']) * TECUns/(LIGHT_SPEED*1e-9)
        if 'C1' in self[prn] and 'C2' in self[prn]:
            return (self[prn]['C2'] - self[prn]['C1']) * TECUns/(LIGHT_SPEED*1e-9)        

    def badness(self, prn):
        '''
        Indicate how `bad' this record is for TEC calculation
        '''

        if prn not in self:
            return True  # Very bad!
        elif 'L1' not in self[prn] or 'L2' not in self[prn]:
            return True
        elif 'C1' not in self[prn] and 'P1' not in self[prn]:
            return True
        elif 'C2' not in self[prn] and 'P2' not in self[prn]:
            return True
        elif abs(self.ctec(prn)) > 990 or abs(self.ptec(prn)) >99999990:
            return True
        elif self[prn]['L2']==0 or self[prn]['P2']==0:
            return True
        return False

class PRNrecord(dict):
    '''
    '''
    
    def __getitem__(self, index):
        '''
        '''
        if index in self:
            return dict.__getitem__(self, index)
        else:
            val = value(0)
            val.lli = ' '
            val.str = ' '
            return val

class records_cnt(object):
    '''
    Class to parse rinex records only for count purpose used with count_rinex
    function.
    '''
    def __init__(self, iscrx):
        self.line = ''
        self.iscrx = iscrx

    def update(self, fid):
        self.line = self.getline(fid)
        self.numrec = toint(self.line[29:32])

    def getline(self, fid):
        if not self.iscrx:
            return fid.next()
        else:
            line = fid.next()
            ln = fid.next()
            try: xyz=line[0]
            except IndexError:
                warn('Duplicate record found, epoch=' + \
                     self.epoch.strftime('%y/%m/%d %H:%M:%S'))
                for n in range(self.numrec+2):
                    line = fid.next()
            if line[0] == '&':
                return line.replace('&', ' ')
            else:
                return ''.join(map(choose, self.line, line))

    def prnlist(self, fid):

        prnlist = []
        if not self.iscrx:
            line = self.line
            for z in range(self.numrec):
                s = z % 12
                if z and not s:
                    line = fid.next()
                prn = btog(line[32 + s * 3]) + '%02d' % \
                                            toint(line[33 + s * 3 : 35 + s * 3])
                prnlist += [prn]
            return prnlist
        else:
            for s in range(self.numrec):
                prn = btog(self.line[32 + s * 3]) + '%02d' % \
                                    toint(self.line[33 + s * 3 : 35 + s * 3])
                prnlist += [prn]
            return prnlist

    def dataline(self, fid, numobs):

        line = fid.next()
        if not self.iscrx and numobs>5:
            line = fid.next()
        return

class records_rnx(object):
    '''
    Parse record headers (epoch lines) in standard RINEX:
    Combine continuation lines if necessary.
    '''
    def __init__(self, baseyear):
        self.line = ''
        self.baseyear = baseyear
        self.epoch = None
        #self.intervals = set()

    def update(self, fid):
        '''
        Process a new epoch line.
        '''
        self.line     = self.getline(fid)
        self.oldepoch = self.epoch
        self.epoch    = recordtime(self.line[:26], self.baseyear)
        self.numrec   = toint(self.line[29:32])
        
        while self.oldepoch and self.epoch and self.epoch<=self.oldepoch:
            warn('Duplicate or earlier record found skiping, epoch= %s' %
                 self.epoch)
            n = 0
            if hasattr(self, 'is_crx'):
                n = 1
            for i in range(self.numrec+n):
                lin = fid.next()
            self.line = self.getline(fid)
            self.epoch = recordtime(self.line[:26], self.baseyear)
            self.numrec   = toint(self.line[29:32])

        #if self.epoch is not None and self.oldepoch is not None:
        #    self.intervals.add(self.epoch - self.oldepoch)
        self.numrec = toint(self.line[29:32])
        self.flag = toint(self.line[28])

    def getline(self, fid):
        return fid.next()

    def prnlist(self, fid):
        '''
        Return the list of PRNs (satellite IDs) included in this epoch line.
        May consume extra lines if there are more than 12 PRNs.
        '''
        prnlist = []
        line = self.line
        for z in range(self.numrec):
            s = z % 12
            if z and not s:
                line = fid.next()
            prn = btog(line[32 + s * 3]) + '%02d' % \
                                        toint(line[33 + s * 3 : 35 + s * 3])
            prnlist += [prn]
        return prnlist

    def get_obs(self, fid, prn, obscodes):
        obsline = PRNrecord()
        for x in range(len(obscodes)):
            ind = x%5
            if not ind:
                line = fid.next()
            val = value(tofloat(line[ind*16 : ind*16+14]))
            val.lli = toint(line[ind*16+14 : ind*16+15])
            val.str = toint(line[ind*16+15 : ind*16+16])
            obsline[obscodes[x]] = val
        return obsline

    def offset(self, fid):
        '''
        Return receiver clock offset (optionally) included at end of epoch line.
        '''
        return tofloat(self.line[68:])

class records_crx(records_rnx):
    '''
    Parse record headers in Compact RINEX:
    each line only contains differences from the previous.
    '''
    def __init__(self, baseyear):
        self.is_crx = True
        self.data = {}
        self.offsetval = None
        records_rnx.__init__(self, baseyear)

    def getline(self, fid):
        self.offsetval = None
        line = fid.next()
        if not line:
            return ''.join(map(choose, self.line, line))
        elif line[0] == '&':
            return line.replace('&', ' ')
        else:
            return ''.join(map(choose, self.line, line))

    def prnlist(self, fid):
        prnlist = []
        for s in range(self.numrec):
            prn = btog(self.line[32 + s * 3]) + '%02d' % \
                                       toint(self.line[33 + s * 3 : 35 + s * 3])
            prnlist += [prn]
        return prnlist

    def get_obs(self, fid, prn, obscodes):
        obsline = {}
        numobs = len(obscodes)

        if prn not in self.data:
            self.data[prn] = {}
            self.data[prn]['arcs'] = [dataArc() for n in range(numobs)]
            self.data[prn]['llis'] = [charArc() for n in range(numobs)]
            self.data[prn]['strs'] = [charArc() for n in range(numobs)]

        line = fid.next()
        vals = line.split(' ', numobs)

        for c, v in enumerate(vals[:numobs]):
            if len(v) >= 2 and v[1] == '&':
                self.data[prn]['arcs'][c] = dataArc(toint(v[0]))
                self.data[prn]['arcs'][c].update(toint(v[2:]))
            elif v.rstrip():
                self.data[prn]['arcs'][c].update(toint(v))
            elif v.rstrip():
                raise ValueError('Uninitialized data arc.')

            obsline[obscodes[c]] = value(self.data[prn]['arcs'][c].get())

        if len(vals) > numobs:
            for c, l in enumerate(vals[numobs][0:numobs*2:2]):
                self.data[prn]['llis'][c].update(l)
                obsline[obscodes[c]].lli = self.data[prn]['llis'][c].get()
            for c, s in enumerate(vals[numobs][1:numobs*2:2]):
                self.data[prn]['strs'][c].update(s)
                obsline[obscodes[c]].str = self.data[prn]['strs'][c].get()
        '''
        for x in range(numobs):
            val = value(self.data[prn]['arcs'][x].get())
            val.lli = self.data[prn]['llis'][x].get()
            val.str = self.data[prn]['strs'][x].get()
            obsline[obscodes[x]] = val
        '''
        return obsline

    def offset(self, fid):
        if self.offsetval is not None:
            return self.offsetval
        line = fid.next()
        if len(line) >= 2 and line[1] == '&':
            self.offsetArc = dataArc(toint(line[0]))
            self.offsetArc.update(toint(line[2:]))
        elif line.rstrip() and 'offsetArc' in self.__dict__:
            self.offsetArc.update(toint(line))
        elif line.rstrip():
            raise ValueError('Uninitialized clock offset data arc.')
        else:
            return 0.
        return self.offsetArc.get()/1000000000

class dataArc(object):
    '''
    Numeric records in Compact RINEX are Nth-order differences
    from previous records.
    Difference order is usually 3. Fields are separated by space.
    LLI and STR are kept separately at the end of the line in one character
    each.
    '''
    def __init__(self, order=3):
        self.order = order
        self.data = []
        self.index = 0

    def update(self, value):
        if self.index < self.order:
            self.data.append(value)
            self.index += 1
        else:
            self.data[self.order - 1] += value

        for diff in range(self.index - 2, -1, -1):
            self.data[diff] += self.data[diff + 1]
        return self.data[0]

    def get(self):
        if len(self.data):
            return self.data[0]/1000.
        else:
            return 0

class charArc(object):
    '''
    LLI and STR records in Compact RINEX only record changes from the previous
    record; space indicated no change.
    '''
    def __init__(self):
        self.data = '0'
    def update(self, char):
        self.data = ''.join(map(choose, self.data, char))
    def get(self):
        return toint(self.data)

class TECData(GPSData):
    '''
    '''
    def __init__(self, filename, secs=None, satsystem='G'):
        header = {}
        station = {}
        self.phasearcs = {}
        self.prns = set()
        self.satsystem = None
        self.rec_bias = 0
        self.bad_prn = []    
        self.fo = open_file(filename)
        station['fullname'] = self.fo.next().strip()
        station['location'] = tuple(self.fo.values(fmt=float)[:3])
        station['receiver'] = self.fo.next().strip()
        header['receivertype'] = station['receiver']
        header['marker'] = {0: filename.split('/')[-1][:4]}
        station['lcode'] = header['marker'][0]
        self.fo.next()
        self.date = GPSDateTime(*self.fo.values(fmt=int))
        ###first record
        prn = '%s%02d' % (satsystem, int(self.fo.values()[1]))
        self.prns.add(prn)
        if secs is None:
            vals = self.fo.values()
            dt = self.date.replace(hour=int(vals[2]), minute=int(vals[3]), second=int(vals[4]))
        else:
            while True:
                vals = self.fo.values()
                dt = self.date.replace(hour=int(vals[2]), minute=int(vals[3]), second=int(vals[4]))
                if dt.second in secs:
                    break
        self.append(TECrecord(dt, prn, vals))
        
        while True:
            try:
                vals = self.fo.values()
            except StopIteration:
                break
            if 'PRN' in vals:
                prn = '%s%02d' % (satsystem, int(vals[1]))
                self.prns.add(prn)
                continue
            self.add(prn, vals, secs)            
        
        self.intervals = set([self[i].epoch-self[i-1].epoch 
                              for i in range(1, len(self))])
        [self.checkbreak(i) for i in range(len(self))]
        self.endphase(len(self)-1, self.prns)
        self.header = header
        self.station = station
        
        
    def add(self, prn, vals, secs=None):
        '''
        '''
        dt = self.date.replace(hour=int(vals[2]), minute=int(vals[3]), second=int(vals[4]))
        if (secs is not None) and (dt.second not in secs):            
            return
        
        if dt<self[0].epoch:
            self.insert(0, TECrecord(dt, prn, vals))
        elif dt>self[-1].epoch:
            self.append(TECrecord(dt, prn, vals))
        else:
            dum = [dt>rec.epoch for rec in self]
            pos = dum.index(False)
            if dt==self[pos].epoch:
                self[pos].add_values(prn, vals)
            else:
                self.insert(pos, TECrecord(dt, prn, vals))
        
    def checkbreak(self, n):
        '''
        '''
        
        for prn in set(self[n]).union(self.phasearcs):
            if n>0:
                if self[n].epoch-self[n-1].epoch>1800:
                    self.breakphase(n, prn)
                    continue
            
            bad = prn not in self[n] 
                                    
            if not bad and prn not in self.phasearcs:
                self.phasearcs[prn] = [[n, None]]
                continue
            if not bad and self.phasearcs[prn][-1][1] is not None:
                self.phasearcs[prn] += [[n, None]]
                continue
            if prn not in self.phasearcs or self.phasearcs[prn][-1][1] is not \
                    None:
                continue
            if bad:
                self.endphase(n, prn)
                continue
            
    def breakphase(self, n, prn):
        '''
        Begin new phase-connected-arc for satellite prn.
        '''
        if isinstance(prn, (list, tuple, set, dict)):
            [self.breakphase(p) for p in prn]
        elif prn not in self.phasearcs:
            self.phasearcs[prn] = [[n, None]]
        elif prn not in self[n]:
            self.endphase(n, prn)
        else:
            self.endphase(n, prn)
            self.phasearcs[prn] += [[n, None]] 
            
    def endphase(self, n, prn):
        '''
        End current phase-connected-arc, if any, for satellite prn.
        Ends arc just before the current record.
        '''

        if isinstance(prn, (list, tuple, set, dict)):
            [self.endphase(n, p) for p in prn]
        elif prn in self.phasearcs and self.phasearcs[prn][-1][1] is None:
            self.phasearcs[prn][-1][1] = n-1          
                
class TECrecord(dict):
    '''
    '''
    def __init__(self, epoch, prn, vals):
        self.epoch = epoch
        self.add_values(prn, vals)
    
    def __getitem__(self, index):
        '''
        Allow you to access GPS satellites, eg record['G13'], as
        simply record[13].  For GLONASS or Galileo, you must use the full code.
        '''
        if index == 'epoch':
            return self.epoch
        else:
            return dict.__getitem__(self, index)
        
    def add_values(self, prn, vals):
        '''
        '''
        vars = ['eqTEC', 'ltime', 'ele', 'lat', 'lon', 'sTEC']
        self[prn] = dict(zip(vars, [float(x) for x in vals[5:]]))

class header(object):
    '''
    For each RINEX header type, this holds a list of field objects
    which are defined in the associated line.
    This is for header values which should only occurr once.
    '''

    class field(object):
        '''
        Describes a value in a RINEX header: variable name, position in the
        line, and how to interpret it.
        '''
        def __init__(self, name, start, stop, convert=mystrip):
            self.name = name
            self.start = start
            self.stop = stop
            self.convert = convert

        def read(self, line):
            return self.convert(line[self.start:self.stop])
        '''
        def __deepcopy__(self, memo={}):
            #Fix deepcopying in Python 2.4.
            newfield = self.__class__(deepcopy(self.name, memo),
                              deepcopy(self.start, memo),
                              deepcopy(self.stop, memo), self.convert)
            memo[id(self)] = newfield
            return newfield
        '''

    def __init__(self, field_args):
        self.values = [header.field(*fargs) for fargs in field_args]

    def read(self, line, header, recordnum, epoch=None):
        for field in self.values:
            header[field.name] = field.read(line)

class listheader(header):
    '''
    For multiple header values

    The value may change for different observation records.
    If multiple instances are at the same record number, the last is used.
    They are accessed by record number; whichever value is valid for that
    record is returned.
    '''
    def read(self, line, header, recordnum, epoch=None):
        for field in self.values:
            if field.name not in header:
                header[field.name] = {}
            if field.name == 'comment':
                if recordnum not in header[field.name]:
                    header[field.name][recordnum] = []
                header[field.name][recordnum] += [field.read(line)]
            else:
                header[field.name][recordnum] = field.read(line)
                if epoch is not None:
                    header[field.name][recordnum].epoch = epoch

RINEX = {
    'CRINEX VERS   / TYPE' : header([('crnxver', 0, 3, crxcheck),
                              ('is_crx', 0, 0, truth)]),
    'CRINEX PROG / DATE  ' : header([('crnxprog', 0, 20),
                              ('crxdate', 40, 60),
                              ('is_crx', 0, 0, truth)]),
    'RINEX VERSION / TYPE' : header([('rnxver', 0, 9, versioncheck),
                              ('filetype', 20, 21, iso),
                              ('satsystem', 40, 41, btog)]),
    'PGM / RUN BY/ DATE  ' : header([('rnxprog', 0, 20),
                              ('run_by', 20, 40),
                              ('filedate', 40, 60)]),
    'PGM / RUN BY / DATE ' : header([('rnxprog', 0, 20),
                              ('run_by', 20, 40),
                              ('filedate', 40, 60)]),
    'COMMENT             ' : listheader([('comment', 0, 60)]),
    'MARKER NAME         ' : listheader([('marker', 0, 60)]),
    'MARKER NUMBER       ' : listheader([('markernum', 0, 20)]),
    'APPROX POSITION XYZ ' : listheader([('markerpos', 0, 42, to3float)]),
    'OBSERVER / AGENCY   ' : header([('observer', 0, 20),
                              ('agency', 20, 60)]),
    'REC # / TYPE / VERS ' : header([('receivernum', 0, 20),
                              ('receivertype', 20, 40),
                              ('receiverver', 40, 60)]),
    'ANT # / TYPE        ' : listheader([('antennanum', 0, 20),
                              ('antennatype', 20, 40)]),
    'ANTENNA: DELTA H/E/N' : listheader([('antennashift', 0, 42, to3float)]),
    'WAVELENGTH FACT L1/2' : listheader([('ambiguity', 0, 53, wavelength())]),
    '# / TYPES OF OBSERV ' : listheader([('obscodes', 0, 60, obscode())]),
    'INTERVAL            ' : header([('interval', 0, 10, tofloat)]),
    'TIME OF FIRST OBS   ' : header([('firsttime', 0, 43, headertime),
                              ('firsttimesys', 51, 54)]),
    'TIME OF LAST OBS    ' : header([('endtime', 0, 43, headertime),
                              ('endtimesys', 51, 54)]),
    'RCV CLOCK OFFS APPL ' :
                      listheader([('receiverclockcorrection', 0, 6, toint)]),
    'LEAP SECONDS        ' : listheader([('leapseconds', 0, 6, toint)]),
    '# OF SATELLITES     ' : header([('numsatellites', 0, 6, toint)]),
    'PRN / # OF OBS      ' :
                      header([('obsnumpersatellite', 3, 60, satnumobs())])
}

def procheader(fid, RINEX, header, numrec, numlines=itertools.repeat(0),
                        epoch=None):
    '''
    Parse Rinex headers into a dictionary.
    '''

    for n in numlines:
        try:
            line = '%-80s' % fid.next()  # pad spaces to 80
        except StopIteration:
            break
        label = line[60:]
        if label.strip() == 'END OF HEADER':
            break
        if label in RINEX:
            RINEX[label].read(line, header, numrec, epoch)
        else:
            warn('Header line ' + label + ' unrecognized; ignoring')

def open_rinex(filename, verbose=False):
    '''
    Open a rinex file (rinex, crinex, gziped y/o tarred), return an utl.open_file
    object
    '''

    if tarfile.is_tarfile(filename):
        if verbose: print 'Unpacking tarfile.'
        zfile = tarfile.open(filename)  # Automatically handles tar.gz,bz2
        zfile = zfile.extractfile(zfile.next())
    elif filename.lower().endswith('.gz'):
        if verbose: print 'Gunzipping file.'
        zfile = gzip.open(filename)
        zfile.name = filename[:filename.rfind('.gz')]
    elif filename.endswith('.Z'):
        zfile = StringIO(Popen(["zcat", filename],
                         stdout=PIPE).communicate()[0])
        #zfile = StringIO(subprocess.check_output(['zcat', filename]))
    else:
        zfile = open(filename)
    if not re.search('\.[0-9]{2}[OoDd]$', filename):
        if verbose: warn(filename + ' Invalid Filename format')
    return open_file(zfile)

def count_rinex(filename):
    '''
    Simplified version of read_rinex that return the number of records of the
    file.
    '''
    fid = open_rinex(filename)
    header = {}
    procheader(fid, RINEX, header, 0)
    record = records_cnt('is_crx' in header)
    count = 0
    while True:
        try:
            record.update(fid)
        except StopIteration:
            break
        count += 1
        for prn in record.prnlist(fid):
            record.dataline(fid, header['obscodes'][0])
    fid.close()
    return count

def read_rinex(filename, check=True, pcode=None, headeronly=False,
               station=False, verbose=False):
    '''
    Parse Rinex file and return a GPSData object with observables data.
    '''
    print 'Parsing file: %s' % filename
    fid = open_rinex(filename, verbose)
    obsdata = GPSData()
    procheader(fid, RINEX, obsdata.header, 0)
    baseyear = obsdata.timesetup(filename)
    if headeronly:
        fid.close()
        return obsdata
    if 'is_crx' in obsdata.header:
        record = records_crx(baseyear)
    else:
        record = records_rnx(baseyear)

    while True:
        try:
            record.update(fid)
        except StopIteration:
            break

        if record.flag == 5:
            procheader(fid, RINEX, obsdata.header, len(obsdata),
                       xrange(record.numrec), record.epoch)
        elif record.flag == 4:
            procheader(fid, RINEX, obsdata.header, len(obsdata),
                       xrange(record.numrec))
        elif 2 <= record.flag <= 3:
            obsdata.inmotion = record.flag == 2
            procheader(fid, RINEX, obsdata.header, len(obsdata),
                       xrange(record.numrec))
        elif 0 <= record.flag <= 1 or record.flag==6:
            obsdata.append(GPSrecord(record, record.offset(fid)))
            for prn in record.prnlist(fid):
                dataline = record.get_obs(fid, prn, obsdata.obscodes())
                obsdata.update(prn, dataline)
            obsdata.checkbreak(record.flag==6)
    fid.close()
    
    if check:
        obsdata.endphase(obsdata.prns)
        obsdata.check()
    if pcode:
        for x in obsdata.header['obscodes']:
            i = obsdata.header['obscodes'][x].index('C1')
            obsdata.header['obscodes'][x][i] = 'P1'

    if station==True:
        station = utl.Stations(False)
        station.read(obsdata.header['marker'][0],
                     xyz2lla(*obsdata.header['markerpos'][0]))
        obsdata.station = station
    elif station<>False:
        obsdata.station = station
    return obsdata

def nvd_to_rnx(input_file, marker='abcd', rnx_path=localpath, date=None,
                obscodes=rnx_codes, interval=1, alt_name=None, **kwargs):
    '''
    Read a Novatel binary file (nvd) and create a Rinex observation file.

    Inputs:
        input_file: binary novatel file
        marker: Four char station name for rinex header.
        rnx_path: path where the rinex file will be created.
        obscodes: list of rinex codes
        interval: interval for decimated data.
        kwargs: aditional information for rinex header.
    Outputs: Rinex filename if created and number of records
    '''

    def decode_id43(buffer):
        '''
        Unpack id43 records from buffer, return  a list of records objects
        '''
        num_obs = unpack('=L', buffer[:4])[0]
        ind = 4
        rec_list = []
        last_rec = ReadRecord()
        for i in range(num_obs):
            record1 = ReadRecord('BB', buffer[ind:ind+2])
            record2 = ReadRecord('=dfdffffL', buffer[ind+4:ind+44])
            record1.vals += record2.vals
            record1.id43()
            #print record1.prn, record1
            if record1.prn <> last_rec.prn:
                last_rec = record1
            else:
                last_rec.update(record1)
                rec_list.append(last_rec)
            ind += 44
        return rec_list

    if input_file.lower().endswith('.nvd.gz') or \
        input_file.lower().endswith('.nvd'):
        try:
            rawfile = open_file(input_file)
            print 'Opening file', input_file
        except IOError:
            print 'Error reading file:', input_file
            return None, 0
    else:
        print input_file, 'Invalid nvd file'
        return None, 0

    nrec = 0
    intervals = set()
    sync = [0, 0, 0]
    last_obs = datetime(1,1,1)

    while True:
        try:
            ch = unpack('B', rawfile.read(1))[0]
        except:
            break
        sync = [sync[1], sync[2], ch]
        #print old_pos, sync
        if sync == [0xAA, 0x44, 0x12]:
            file_pos = rawfile.tell()
            buff0 = rawfile.read(1)
            hdr_len = unpack('B', buff0)[0]
            buff1 = rawfile.read(hdr_len-4)
            header = ReadRecord('=HBBHHBBHlLHH', buff1)
            header.id43_header()
            week_sec = header.milli_sec/1000
            gps_sec = header.week*SECONDS_IN_WEEK + week_sec
            dow = week_sec/86400
            daysec = week_sec%86400
            hh = daysec/3600
            mm = (daysec%3600)/60
            ss = daysec%60
            rec_date = gdate(header.week*7 + dow + jday_gps0) + (hh, mm, ss)
            buff2 = rawfile.read(header.msg_len)
            buffer = '\xaa\x44\x12'+buff0+buff1+buff2
            crc_1 = unpack('=L', rawfile.read(4))[0]
            crc_2 = CRC32(buffer) & 0xFFFFFFFF

            if crc_1 <> crc_2:
                rawfile.seek(file_pos, 0)
                print 'CRC_error', crc_1, crc_2, rec_date
                continue

            is_dec = gps_sec%interval == 0

            if header.msg_id == 43 and is_dec:
                obs_list = decode_id43(buff2)
                if nrec == 0 and len(obs_list) <> 0:
                    first_obs = GPSDateTime(*rec_date)
                    date_obs  = GPSDateTime(*rec_date)
                    if alt_name:
                        rnx_name = '%s.%02do' % (alt_name, date_obs.year%100)
                    else:
                        rnx_name = '%s%03d0.%02do' % (marker, date_obs.doy,
                                                      date_obs.year%100)
                    print 'GPS week set to:',header.week
                    print 'Creating rinex file:', os.path.join(rnx_path, rnx_name)
                    print 'Given obscodes:', obscodes
                    print 'Available obscodes:', obs_list[0].keys()
                    obscodes = [x for x in obscodes if x in obs_list[0]]
                    Rinex = MkRinex(marker, interval, obscodes)
                    file_tmp = os.tmpfile()
                elif nrec <> 0 and len(obs_list) <> 0:
                    intervals.add((datetime(*rec_date)-date_obs).seconds)
                    date_obs = datetime(*rec_date)
                    if date_obs<=last_obs:
                        print 'Earlier data found'
                        continue
                    last_obs = date_obs
                else:
                    continue
                prn_list = ['G%02d' % rec.prn for rec in obs_list]
                sat_line = Rinex.mk_satline(date_obs, prn_list)
                dat_line = ''
                for rec in obs_list:
                    obs_values = [rec[x] for x in obscodes]
                    lli_values = [0 for x in obs_values]
                    snr_values = [Rinex.snr_flag(rec, x) for x in obscodes]
                    dat_line += Rinex.mk_recordline(obs_values, lli_values,
                                snr_values)
                file_tmp.write(sat_line+dat_line)
                nrec += 1
    if nrec<>0:
        kwargs['endtime'] = last_obs
        if intervals:
            if min(intervals) <> interval:
                print 'Correcting interval %d secs. to %d secs.' % \
                    (interval, min(intervals))
                Rinex.interval = min(intervals)
        header = Rinex.mk_header(first_obs, **kwargs)
        if not os.path.exists(rnx_path):
            os.makedirs(rnx_path)
        file_tmp.seek(0,0)
        file_rnx = open(os.path.join(rnx_path, rnx_name), 'w')
        file_rnx.write(header)
        file_rnx.write(file_tmp.read())
        rawfile.close()
        file_tmp.close()
        file_rnx.close()
        print 'Rinex records:', nrec
        return os.path.join(rnx_path, rnx_name), nrec
    else:
        print 'No records found for given date'
        return None, 0

def lb2_to_rnx(input_file, marker='site', rnx_path=localpath, date=None,
                obscodes=rnx_codes, interval=1, alt_name=None, wpos=False,
                **kwargs):
    '''
    Read a Leica file lb2 and create a Rinex observation file.

    Inputs:
        input_file: binary leica file (id37 or lb2)
        marker: Four char station name for rinex header.
        rnx_path: path where the rinex file will be created.
        obscodes: list of rinex codes
        interval: interval for decimated data
        kwargs: aditional information for rinex header.
    Outputs: Rinex filename if created and number of records
    '''

    def decode_id3(buffer, ndata):
        '''
        Unpack id3 records from buffer, return  a list of records objects
        '''
        ind = 6
        obs_list = []
        while ind<ndata:
            if len(buffer[ind:ind+3])<>3:
                ind += len(buffer[ind:ind+3])
                continue
            record = ReadRecord('=Bh', buffer[ind:ind+3])
            record.id3(buffer, ind+3)
            if sum(record.values())<>0:
                obs_list.append(record)
            ind = record.ind
            #print record.ind
        if ind==ndata:
            return obs_list
        else:
            return None

    def decode_id4(buffer, ndata):
        '''
        Unpack id4 records from buffer, return  a list of records objects
        '''
        nsat = unpack('B', buffer[0])[0]
        ind = 1
        pos_list = []
        while ind<ndata:
            record = ReadRecord('=Bhh', buffer[ind:ind+5])
            record.id4()
            pos_list.append(record)
            ind += 5
        if ind==ndata:
            return pos_list
        else:
            return None

    def decode_id85(buffer, ndata):
        '''
        Unpack id85 records from buffer, return  a list of records objects
        '''
        record = ReadRecord('B', buffer[0])
        record.id85(buffer)
        #if record.ind == ndata:
        return record

    def decode_id37(buffer, header, ind):
        '''
        Unpack id37 records from buffer, return  a list of records objects
        '''
        obs_list = []
        for i in range(header.nprn):
            record = ReadRecord('=Bdhdh', buffer[ind:ind+21])
            record.id37()
            if sum(record.values())<>0:
                obs_list.append(record)
            ind += 24
        return obs_list

    if input_file.lower().endswith('.lb2') or \
        input_file.lower().endswith('.lb2.gz'):
        try:
            rawfile = open_file(input_file)
            print 'Opening file', input_file
        except IOError:
            print 'Error reading file:', input_file
            return None, 0
    else:
        print input_file, 'Invalid lb2 file'
        return None, 0

    file_name = input_file.split('/')[-1]
    if date:
        file_date = GPSDateTime(date)
    else:
        try:
            file_date = GPSDateTime(file_name.split('.lb2')[0])
        except:
            print 'Filename does not have a valid date stamp: yymmdd'
            return None, 0

    week = file_date.GPSweek
    print 'GPS week set to:', week
    nrec = 0
    prec = 0
    sync = [0, 0]
    intervals = set()
    first_jd = 0
    first_obs = 0
    last_obs = datetime(1,1,1)
    
    while True:
        try:
            ch = unpack('B', rawfile.read(1))[0]
        except:
            break
        sync = [sync[1], ch]
        if sync == [0x9C, 0xAE]:
            file_pos = rawfile.tell()
            buff0 = rawfile.read(2)
            rec_len = unpack('H', buff0)[0]
            if rec_len >8192:
                rawfile.seek(file_pos, 0)
                continue
            buff1 = rawfile.read(1)
            rec_id = unpack('B', buff1)[0]
            buff2 = rawfile.read(rec_len-5)
            buff3 = rawfile.read(2)
            checksum1 = unpack('h', buff3)[0]
            num_bytes = 3+rec_len-5
            checksum2 = sum(unpack('%iB' % num_bytes, buff0+buff1+buff2))*-1

            '''
            if checksum1 <> checksum2:
                print 'Checsum error'
                rawfile.seek(file_pos, 0)
                continue
            '''

            save_rnx = False
            save_pos = False
            ndata = rec_len-5

            if rec_id == 0x02:
                #obs_list = decode_id2(buff2, ndata)
                continue
            elif rec_id == 0x03:
                if ndata <=6:
                    rawfile.seek(file_pos, 0)
                    continue
                header = ReadRecord('=3chB', buff2[:6])
                header.id3_header()
                is_dec = header.gpstow%(interval*10) == 0
                week_sec = header.gpstow/10
                dow = week_sec/86400
                day_sec = week_sec % 86400
                hms = (day_sec/3600, (day_sec % 3600)/60, day_sec % 60)
                rec_date = gdate(week*7 + dow + jday_gps0) + hms

                obs_list = decode_id3(buff2, ndata)
                save_rnx = True
                save_pos = False
            elif rec_id == 0x04:
                pos_list = decode_id4(buff2,ndata)
                save_pos = True
                save_rnx = False
            elif rec_id == 0x85:
                record = decode_id85(buff2, ndata)
                if record.date[2]>80:
                    yr = record.date[2]+1900
                else:
                    yr = record.date[2]+2000
                date_id85 = (yr, record.date[1], record.date[0])
                jd_id85 = jd(*date_id85)
                continue
            elif rec_id == 0x37:
                station_prefix = buff2[:2]
                i0 = 2
                if buff2[3:6] == 'ver':
                    print 'ver True'
                    i0 += 10
                header = ReadRecord('11B', buff2[i0:i0+11])
                header.id37_header()
                is_dec = header.gps_sec%interval == 0
                rec_date = header.gps_date+header.dhms[1:]

                obs_list = decode_id37(buff2, header, i0+13)
                save_rnx = True
                save_pos = False
            else:
                #print 'Lb2 id %d unknown, skipping' % rec_id
                continue

            if save_rnx and is_dec:
                if nrec == 0 and len(obs_list)<>0:
                    first_obs = GPSDateTime(*rec_date)
                    date_obs  = GPSDateTime(*rec_date)
                    if date_obs<file_date:
                        continue
                    elif date_obs>(file_date+timedelta(2)):
                        continue
                    if alt_name:
                        rnx_name = '%s.%02do' % (alt_name, date_obs.year%100)
                    else:
                        rnx_name = '%s%03d0.%02do' % (marker, date_obs.doy,
                                                      date_obs.year%100)
                    print 'Creating rinex file:', os.path.join(rnx_path,
                                                                rnx_name)
                    print 'Given obscodes:', obscodes
                    print 'Available obscodes:', obs_list[0].keys()
                    obscodes = [x for x in obscodes if x in obs_list[0]]
                    Rinex = MkRinex(marker, interval, obscodes)
                    file_tmp = os.tmpfile()
                elif nrec <> 0 and len(obs_list)<>0:
                    intervals.add((datetime(*rec_date)-date_obs).seconds)
                    date_obs = datetime(*rec_date)
                    if date_obs<=last_obs:
                        print 'Earlier data found'
                        continue
                    last_obs = date_obs
                else:
                    continue

                prn_list = ['G%02d' % rec.prn for rec in obs_list]
                sat_line = Rinex.mk_satline(date_obs, prn_list)
                dat_line = ''
                for rec in obs_list:
                    obs_values = [rec[x] for x in obscodes]
                    lli_values = [0 for x in obs_values]
                    snr_values = [Rinex.snr_flag(rec, x) for x in obscodes]
                    dat_line += Rinex.mk_recordline(obs_values, lli_values,
                                snr_values)
                file_tmp.write(sat_line+dat_line)
                nrec += 1

            if save_pos and wpos:
                if prec==0:
                    pos_name = '%s_%02d%02d%02d.pos' % (marker, rec_date[0]%100,
                        rec_date[1], rec_date[2])
                    pos_path = os.path.join(os.path.dirname(rnx_path), 'posit')
                    print 'Creating position file:', os.path.join(pos_path,
                                                                    pos_name)
                    if not os.path.exists(pos_path):
                        os.makedirs(pos_path)
                    pos_file = open(os.path.join(pos_path, pos_name), 'w')
                    pos_file.write('HH MM SS PRN  az el')
                pos_line = ''
                for rec in pos_list:
                    pos_line += '%02d %02d %02d %3d %3d %2d\n' % (rec_date[3],
                                rec_date[4], rec_date[5],rec.prn, rec['az'],
                                rec['el'])
                pos_file.write(pos_line)
                prec +=1
    if prec<>0:
        print 'pos records:', prec
        pos_file.close()

    if nrec<>0:
        kwargs['endtime'] = last_obs
        if intervals:
            if min(intervals) <> interval:
                print 'Correcting interval %d secs. to %d secs.' % \
                    (interval, min(intervals))
                Rinex.interval = min(intervals)
        header = Rinex.mk_header(first_obs, **kwargs)
        if not os.path.exists(rnx_path):
            os.makedirs(rnx_path)
        file_tmp.seek(0,0)
        file_rnx = open(os.path.join(rnx_path, rnx_name), 'w')
        file_rnx.write(header)
        file_rnx.write(file_tmp.read())
        rawfile.close()
        file_tmp.close()
        file_rnx.close()
        print 'Rinex records:', nrec
        return os.path.join(rnx_path, rnx_name), nrec
    else:
        print 'No records found for given date'
        return None, 0

def obs_to_rnx(input_file, marker='site', rnx_path=localpath, date=None,
                obscodes=rnx_codes, interval=1, alt_name=None, pcode='P1', **kwargs):
    '''
    Read a observable ascii file (obs) and create a Rinex observation file.

    Inputs:
        input_file: observables file (GPS-Scinda)
        marker: Four char station name for rinex header.
        rnx_path: path where the rinex file will be created.
        obscodes: list of rinex codes
        interval: interval for decimated data.
        kwargs: aditional information for rinex header.
    Outputs: Rinex filename if created and number of records
    '''

    def decode_obs(rawfile, last, pcode):
        '''
        Unpack obs records from rawfile, return  a list of records objects
        '''
        first = False
        obs_list = last[-1:]
        if len(obs_list)<>0:
            millisec = obs_list[-1].millisec
        while True:
            buff = rawfile.readline()
            if not buff:
                last = []
                break
            dataline = ReadRecord(None, buff)
            if len(dataline.vals)<>11: continue
            dataline.obs_line(pcode)
            if len(obs_list)<>0:
                if millisec < dataline.millisec:
                    last = [dataline]
                    break
            millisec = dataline.millisec
            if dataline.prn not in [dl.prn for dl in obs_list]:
                obs_list.append(dataline)
        return obs_list, last

    if input_file.lower().endswith('.obs') or \
        input_file.lower().endswith('.obs.gz'):
        try:
            rawfile = open_file(input_file)
            print 'Opening file', input_file
        except IOError:
            print 'Error reading file:', input_file
            return None, 0
    else:
        print input_file, 'Invalid obs file'
        return None, 0

    file_name = input_file.split('/')[-1]
    if date:
        file_date = date
    else:
        try:
            file_date = GPSDateTime(file_name)
        except:
            print 'Filename does not have a valid date stamp: yymmdd'
            return None, 0

    week = file_date.GPSweek
    print 'GPS week set to:',week, file_date
    nrec = 0
    millisec = 0
    intervals = set()
    last = []
    last_obs = datetime(1,1,1)

    while True:
        obs_list, last = decode_obs(rawfile, last, pcode)
        #print obs_list, last
        if len(last)==0:
            break
        week_sec = obs_list[0].millisec/1000
        gps_sec = week*SECONDS_IN_WEEK + week_sec
        dow = week_sec/86400
        daysec = week_sec%86400
        hh = daysec/3600
        mm = (daysec%3600)/60
        ss = daysec%60
        rec_date = gdate(week*7 + dow + jday_gps0) + (hh, mm, ss)
        if not gps_sec%interval == 0:
            continue

        if nrec == 0 and len(obs_list) <> 0:
            first_obs = GPSDateTime(*rec_date)
            date_obs  = GPSDateTime(*rec_date)
            if date_obs<file_date:
                continue
            elif date_obs>(file_date+timedelta(1)):
                continue
            if alt_name:
                rnx_name = '%s.%02do' % (alt_name, date_obs.year%100)
            else:
                rnx_name = '%s%03d0.%02do' % (marker, date_obs.doy,
                                              date_obs.year%100)
            print 'Creating rinex file:', os.path.join(rnx_path, rnx_name)
            print 'Given obscodes:', obscodes
            print 'Available obscodes:', obs_list[0].keys()
            obscodes = [x for x in obscodes if x in obs_list[0]]
            Rinex = MkRinex(marker, interval, obscodes)
            file_tmp = os.tmpfile()
        elif nrec <> 0 and len(obs_list) <> 0:
            intervals.add((datetime(*rec_date)-date_obs).seconds)
            date_obs = datetime(*rec_date)
            if date_obs<=last_obs:
                print 'Earlier data found'
                continue
            last_obs = date_obs
        else:
            continue

        prn_list = ['G%02d' % rec.prn for rec in obs_list]
        sat_line = Rinex.mk_satline(date_obs, prn_list)
        dat_line = ''
        for rec in obs_list:
            obs_values = [rec[x] for x in obscodes]
            lli_values = [0 for x in obs_values]
            snr_values = [Rinex.snr_flag(rec, x) for x in obscodes]
            dat_line += Rinex.mk_recordline(obs_values, lli_values, snr_values)
        file_tmp.write(sat_line+dat_line)
        nrec += 1
    if nrec<>0:
        kwargs['endtime'] = last_obs
        if intervals:
            if min(intervals)  <> interval:
                print 'Changing interval', interval, \
                        'to min interval bettwen epochs', min(intervals)
                Rinex.interval = min(intervals)
        header = Rinex.mk_header(first_obs, **kwargs)
        if not os.path.exists(rnx_path):
            os.makedirs(rnx_path)
        file_tmp.seek(0,0)
        file_rnx = open(os.path.join(rnx_path, rnx_name), 'w')
        file_rnx.write(header)
        file_rnx.write(file_tmp.read())
        rawfile.close()
        file_tmp.close()
        file_rnx.close()
        print 'Rinex records:', nrec
        return os.path.join(rnx_path, rnx_name), nrec
    else:
        print 'No records found for given date'
        return None, 0

if __name__=='__main__':
    #print timedelta(1./24*1./2).seconds
    
    rnx1 = read_rinex('/data/users/gps/Downloads/laya2740.09d.Z')#, station=True)
    #data = TECData('/data/users/gps/Downloads/laf2_120418.dat.gz',[0])
    #print rnx1
    #rnx2 = read_rinex('/home/jespinoza/Downloads/laya2101.12o')#, station=True)
    #print rnx2
    #rnx1.merge(rnx2)
    print "done"
    #print data.date
    #data.plot_tec()
    rnx1.calctec(path='/data/users/gps/')
    rnx1.plot_tec()
